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Abstract 

Transonic viscous-inviscid interaction is 
considered using the Euler and inverse compress- 
ible turbulent boundary-layer equations. Certain 
improvements in the inverse boundary-layer method 
are mentioned, along with experiences in using 
various Runge-Kutta schemes to solve the Euler 
equations. Numerical conditions imposed on the 
Euler equations at a surface for viscous-inviscid 
interaction using the method of equivalent sources 
are developed, and numerical solutions are pre- 
sented and compared with experimental data to 
illustrate essential points. 

I. Introduction 

Viscous-inviscid interaction is an important 
and difficult problem in transonic aerodynamics. 
Unfortunately, numerical solutions of the Navier- 
Stokes equations are not presently a practical 
method for routinely solving such problems due to 
computer resource requirements. Consequently, 
much research has been done and must is still 
going on with regard to coupling inviscid and 
viscous flow solvers for treating viscous-inviscid 
interaction. Lock^^^ and Melnik^^) have reviewed 
interaction methods. For the most part, these 
methods consist of using potential flow inviscid 
solution methods and attached flow viscous solu- 
tion methods. Inverse boundary-layer methods are 
being used in some instances (see Le Balleur(^) 
for a review) in order to include separated flow. 

Computational fluid dynamics has recently 
matured to the point that numerical solution of 
the Euler equations can be considered for solving 
two- and three-dimensional flow problems . (^”6) 
cause the Euler equations can handle rotational 
flow, these equations offer more information and 
an extended Mach number range compared to the 
potential flow equations. There has, as yet, not 
been a grea t deal of effort devoted to coupling 
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the Euler equations with viscous flow solvers. The 
work that has been done includes Refs. 7-9, where 
the Euler equations were coupled with a compress- 
ible turbulent inverse integral boundary-layer 
method in order to handle rotational flow that 
may Contain regions of separated flow. The pur- 
pose of this paper is to present further results 
of work involving the Euler and inverse boundary- 
layer equations. These results include: (1) im- 

provements in the inverse boundary-layer method, 
(2) numerical .experiments with regard to Euler 
equation boundary conditions, (3) experience using 
second-order Runge-Kutta schemes with various num- 
ber of stages to solve the Euler equations, (4) 
numerical conditions imposed on the Euler equa- 
tions at a surface and in a wake for viscous- 
inviscid interaction using the equivalent source 
method, (5)displacement surface versus the equiv- 
alent source method of interaction, and (6) numer- 
ical and experimental comparisons. 

II. Viscous Method 

The viscous flow solution method is an inverse 
(meaning the pressui^e distribution is obtained as 
part of the solution rather than being specified as 
in a direct method) integral compressible turbu- 
lent boundary-layer method. This inverse method 
is an extension of the direct method described in 
Ref. 11. Both methods solve the momentum and 
mean-flow kinetic energy integral equations. A 
fourth-order four-stage explicit Runge-Kutta 
scheme is used to solve the inverse equations. 

A distinguishing feature of the direct and in- 
verse integral methods in Refs. 11 and 7 was that 
the dissispation integral 

''o w ^ 

was numerically evaluated at each streamwise loca- 
tion as opposed to using an empirical dissipation 
relation. This was accomplished by using a con- 
stant laminar plus turbulent shear stress in the 
region just at the wall, a Cebeci-Smith type model 
in the inner and outer regions, and the derivative 
of the velocity profile expression valid for 
0 < y < T2) Although this placed a stringent 

requirement on the accuracy of the velocity pro- 
file expression, the method gave good results 
even better than finite difference methods for 
transonic flow over adiabatic surf aces. 
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However, the numerical evaluation of the dissipa- 
tion integral at each streamwise location made 
this integral method relatively slow (with regard 
to computational time) as compared to other inte- 
gral methods. The computational time was not a 
severe limitation for steady two-dimensional flow. 
However, with the extension of this method to un- 
steady two- and three-dimensional flow, it was de- 
sirable to eliminate the need for numerically 
evaluating Eq. (1) at each grid point. In this 
connection, Donegan (^ 0 ) succeeded in correlating 
D, as given by Eq. (1), in terms of the local edge 
Mach number, shape factor, and skin friction co- 
efficient (or shape factor and Reynolds number 
based on momentum thickness). Recently, Thomas^^^^ 
has made improvements in the turbulence model used 
in Eq. (1), particularly near the separation point, 
and Donegan and Thomas have improved the correla- 
tion for D given in Ref. 10. The result of using 
an analytical correlation as opposed to numerical- 
ly evaluating Eq. (1) is an increase in speed of 
0 ( 10 ). 

III. Inviscid Method 

Finite volume spatial discretization is ap- 
plied to the integral form of the time-dependent 
Euler equations and the resulting equations were 
solved using second-order Runge-Kutta time-step- 
ping schemes with various number of stages. Dis- 
sipative terms composed of a blend of second and 
fourth differences are used in this central diff- 
erence scheme and these terms are held constant 
during each stage of the Runge-Kutta solution . 
Convergence to a steady state is accelerated by 
the addition of a forcing term that depends on the 
difference between the local total enthalpy and 
the freestream value of enthalpy. Convergence is 
also accelerated by using a local time step deter- 
mined by the maximum Courant number. Far field 
boundary conditions are based on a characteristic 
combination of variables, and pressure at the wall 
is determined using the normal momentum relation. 
With the exception of the use of second-order 
Runge-Kutta schemes with various number of stages 
and frozen dissipation, the numerical method is 
that of Jameson, Schmidt, and Turkel.(^) 

An advantage of this type of explicit scheme 
is that stability can be achieved for Courant num- 
bers greater than one. By using different stage 
Runge-Kutta schemes, the stability region can be 
expanded (see Fig. 1) and the maximum attainable 
Courant number can be increased as shown in Table 
1. Although a larger Courant number can be 
achieved by an increase in the number of stages, 
the increase in work associated with the increase 
in stages eventually reaches a point of diminish- 
ing returns. For example, a scheme with a small 
value of R (see Table 1) should probably be used 
in the early cycles and a scheme with a large val- 
ue of R thereafter. Numerical experiments indi- 
cate that the four-stage scheme is a reasonable 
compromise. The use of second-order accurate 
scheme in time as compared to a fourth-order 
scheme as used in Ref. 6 has the advantage of re- 
quiring slightly less storage. Also, because 
steady state solutions are of interest here, and 
because no noticeable improvement was found in the 
results using a fourth-order scheme as compared to 
a second-order scheme, the method used was the 
second-order four-stag^e scheme with a maximum 
Courant number of 2/2 . 
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Fig. 1 Stability Region for Various Stage 
Second-Order Runge-Kutta Schemes 

Table 1. Second-Order R-Stage Runge-Kutta 
(minimal storage) 
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.76 

.63 1/5, 1/5, 1/3, 1/2 

6 
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.64 1/7, 1/7, 1/4, 1/3, 1/2 
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efficiency for zero dissipation 
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- CFL 
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efficiency for frozen dissipation 


IV. Viscous-Inviscid Coupling 


The displacement surface concept where the 
inviscid solution is carried out on a grid that 
is displaced from the actual body by the amount 
of the boundary-layer displacement thickness, 

6 , is the most commonly used method of viscous- 
inviscid interacton. This approach, however, 
requires that a new grid be generated after each 
boundary-layer solution. A viscous-inviscid 
interaction approach that does not require a new 
grid to be generated after each boundary-layer 
solution is the method of equivalent sources of 
Lighthill . In this method, information from 
the viscous solution is used to specify a distri- 
bution of sources (either positive or negative) 
on the surface and in the wake, and this source 
distribution is used as a boundary condition in 
the inviscid solution. Assuming no attempt is 
made to align some portion of the grid with the 
wake, only one grid must be generated. Unlike a 
potential flow and boundary-layer interaction 
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method, an Euler equation and boundary-layer 
interaction method that uses the equivalent source 
concept requires that information be specified 
for the additional equations of momentum and 
energy. Development of the information necessary 
to use the equivalent source concept with an Euler 
equation and boundary-layer interaction method 
follows* This develonment is based on the work of 
Johnston and Sockol.^^^^ Their work is reviewed 
and then specific relations for the elements of. 
the g vector of the Euler equations at a surface 
are obtained. 

To illustrate the approach consider the 
steady two-dimensional Navier-Stokes equations in 
cartesian coordinates x, y 


if . 

8x 3y 


0 


and the steady two-dimensional Euler equations 


+ if 

3x 3y 


0 


(3) 


— - 
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pv 

puu + p 


puv 


g =• 


puv 


pw + p 

(e + p)u 


(e + p)v 
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_ _ 


and T is the shear stress and q is the heat flux. 
Using the composite function for F, and a similar 
one for G (a point not mentioned in Ref. 16); 

Eq. (A) becomes 



Using Eq. (6) and the definitions of g, f, and 
g, the following conditions on the elements of g^ 
are obtained. 

The term (pv)^ is given by 
a fh 

(pv) = (pv) + -r- [(pu) - pujdy (7) 

O O dX ( O 

' o 

For no porosity in the boundary-layer solution 
((pv)^ = 0] 

( 8 ) 


(9) 
' o 

The term (puv) is given by 
o 

(puv) = (puv - t) + 
o o 

[(pu^ + p) - (pu^ + p)]dy (10) 

o 
o 


3x 


(pv) 




where 6 is defined as 


(P 


* r 

u) 6 « [(pu) - pu]dy 

a o 



and u, V are velocity components in the x, y 
directions, and p, p, and e are the pressure, 
density and total energy per unit volume. An 
explicit description of the elements of F and ^ 
are not needed. Integrating Eqs. (2) and (3) 
with respect to y over 0 < y < h, and considering 
the solution vectors g and 5 to coincide for 
y > h (where h is taken outside the viscous 
region)- the two integrals can be combined to 

obtain(i8) 

o 

where the subscript o indicates y * 0. To avoid 
solving the Navier-Stokes equations, the exact 
solution F is represented^^°^ by.^ composite ^ 
function F , where F^-F=f+f-f , and f is 
a solution^of the boundary-layer equa?ions 

where 


“ “ 

• 

pv 

pu 


puu + p 


puv - T 


i = 
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P 

(e + p)u 


(e + p)v - UT - q 



_ _ 


For no-slip boundary conditions for the boundary- 
layer solution (u =* 0) , and taking the boundary- 
layer pressure equal to the pressure from the 
Euler solution at the surface 

(puv) = * T + ^ [(pu^) (6 + 0)] (11) 

o o ax o 

where 0 is defined as 

(pu^) (6 + 0) = [ I(pu^) - pu^]dy (12) 

‘'o 

As pointed out in Ref. 16 this approach will 
not provide the information necessary to obtain 
the pressure, and a specific approach to obtain 
the third element of g is not given in Ref. 16. 
The pressure is obtained here through an extension 
of the work of Rizzi(^^^ by* including a surface 
porosity term in Rizzi’s normal momentum relation. 
This relation is derived by Thomas^^^^ and the 
influence of including or neglecting the porosity 
term is demonstrated in the next section. The 
term (pv*^ + p)^, therefore, is obtained by deter- 
mining p as mentioned, and determining (pv ) by 
Eq. (8)wSere the density is obtained from the 
previous time step. 

The term ( (e + p)v]^ is given by 

[(e + P)v]^ * [(e + p)v - UT - q]^ + 

3 

^ {[(e + P)u]^ - [(e + p)u]}dy (13) 

-'o 

Using no-slip and no porosity boundary_cond it ions 
for the boundary-layer solution (u^ “ * 
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an adiabatic surface [ (q)^ * 0], and the defini- 
tion of total enthalpy (pH e + p) , Eq. (13) be- 

comes 

d fh 

I(e + p)v]^ - (pvH)^ = — I ((puH)^ - puH]dy 

' o 

(14) 

The boundary-layer method was developed for 
an adiabatic surface with variable total enthalpy 
across the boundary layer that takes into account 
total enthalpy overshoot and nonunity Prandtl 
number. (18) A correlation for the integral in 
Eq. (14) has not been developed as yet, hence the 
approximation H is taken to prevent having to 

numerically evaluate Eq. (14) at each point. This 
approximation yields 

I(e + p)v]^ = (pvH)^ “ Ho ^ I(PU)^«*] 


(15) 

which, by Eq. (8), is now simply an identity. 

It is interesting to note that the combina- 
tion of Eqs. (8) and (11) produce the von Karman 
momentum integral equation. Hence, the results 
of this section can be summarized as 


(pv) = ^ 
o dx 

[(pu)^ 6 ] 

(8) 

(puv) = u 
o o 


(16) 

V 2. 1 

(pv ) = — 

[( pu )^ 6*)}2 

(17) 

o 



[(e + p)v] 

o 

= (pvH) = H ^ [(pu) 6*] 
o o dx o 

(15) 


and the pressure, p , is determined by the above 
mentioned extenison of Rizzi’s method, (^^) 


V. Results 

Numerical results are compared in Refs, 8 
and 19 with experimental data taken on the RAE 
2822 airfoil(20) that show good agreement. ' The 
computations of Refs. 8 and 19 were carried out 
at the geometric angle of attack, cu , of the ex- 
periment as opposed to the corrected angle of 
attack, a^, suggested in Ref. 20 to account for 
wall interference. Recent numerical experiments 
conducted to investigate the sensitivity of the 
solution to the grid, indicate a rather surprising 
sensitivity of lift to the location of the far 
field boundary as indicated in Fig. 2. The 
results in Fig. 2 were obtained by changing the 
location of the far field boundary, while main- 
taining the same far field boundary condit ions(8) 
until there was no further change in the solution. 
Using the grid with the far field boundary located 
such that no change in the. solution due to the 
grid would be expected, the computations for the 
RAE 2822 airfoil for = 0.730, n * 3.19**, and 
Re (freestream Reynolds number basld in chord) = 
6.5 X 10^ were repeated. These results are pre- 
sented in Fig. 3 for both the geometric angle of 
attach (a = Og = 3.19") and the corrected angle 
of attack (a - = 2.78®) suggested by the 

experimenters. (20) as can be seen in Fig. 3, the 
agreement between the computations and the experi- 
mental data for a = a = 2.78" is better than for 
c 


Q “ ttg = 3.19". The freestream Mach number cor- 
rection of 0.004 used in Fig. 3 for the numerical 
solutions was that used by Lock. It appears, 
therefore, that in view of the good agreement be- 
tween numerical and experimental results obtained 
in Fig. 3 by accounting for the sensitivity of the 
far field boundary and using the corrected angle 
of attack, the good agreement obtained previous- 
ly(8, 19)' using a = Og = 3.19" was fortuitous. 
Further experimental results without wall inter- 
ference, or with minimal wall interference and 
accurate far field measurements are needed. 
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Fig. 2 Influence of the Far Field Boundary 
Location on the Lift Coefficient for 
the RAE 2822 Airfoil at * 0.734 
and a = 3.19" (Inviscid) 




Fig. 3 Viscous-Inviscid Interaction Results 
for the RAE 2822 Airfoil Using Geo- 
metric (3.19") and Exper imentor *s(^^) 
Suggested Corrected (2.78") Angle of 
Attack 
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The normal momentum relation derived by 
to obtain surface pressure was based 
on an impermeable surface. Thomas has extend- 
ed this work to include a permeable surface for 
viscous-inviscid interaction. The numerical re- 
sults in Fig. 3 included this new normal momentum 
relation with a permeable surface. A comparison 
of numerical results obtained with and without 
the permeable surface term is given in Fig. A for 
the same flow conditions as Fig. 3. The results 
in Fig. 4 indicate that the influence of this 
term is small, although the influence the term 
does have is to improve the agreement with experi- 
ment slightly on the upper surface at the begin- 
ning of the shock and in the aft region of the 
lower surface. 




Fig. 4 Viscous-inviscid Interaction Results 
for the RAE 2822 Airfoil With and 
Without the Porosity Term in the Sur- 
face Pressure Boundary Condition 


Because the displacement surface method of 
viscous-inviscid interaction is the most commonly 
used method of calculation, a comparison is pre- 
sented in Fig. 5 of the displacement surface 
method and the method of equivalent sources. The 
flow conditions used to obtain the results in Fig 
5 are the same as used to obtain the results in 
Figs- 3 and 4. Also, the normal momentum rela- 
tion allowing for a permeable surface was used. 
There is some difference between the two methods 
of performing interaction computations as in- 
dicated in Fig- 5. The difference in shock 
location, for example, is of the order of the 
distance between grid points in this region. 

As mentioned, the method of equivalent sources 
requires that only one grid be generated, where- 
as, the displacement surface requires a new grid 
for each new boundary-layer displacement surface. 
The method of equivalent sources has been found 
the easiest to use once all the source relations 
are derived and coded. 



Fig. 5 Viscous-inviscid Interaction Results 
for the RAE 2822 Airfoil Using Dis- 
placement Surface and Equivalent Source 
Methods of Interaction 


VI . Conclud ing Remarks 


The results presented involved improvements 
and experiences with previous work in using Euler 
and inverse boundary-layer equations for treating 
transonic viscous-inviscid interaction. Improve- 
ments in the inverse boundary-layer method includ- 
ed the handling of the turbulence modeling, par- 
ticularly near separation, and a correlation for 
the dissipation integral which eliminated the 
need for numerical integration, and thereby re- 
duced the computational time of the viscous sol- 
utions. Solutions of the Euler equations indi- 
cated a significant influence of the far field 
boundary on the lift of a supercritical airfoil 
with a reasonably strong shock on the upper sur- 
face. The computed lift did not change once the 
far field boundary was moved far from the airfoil. 
This observation is receiving further investiga- 
tion. The use of second-order Runge-Kutta schemes 
with various number of stages indicated that a 
second-order four-stage scheme with a maximum 
Courant number of 2/2 was a reasonable compromise 
for solving the Euler equations. Additional num- 
erical surface conditions for the momentum and 
energy equations were developed and used in the 
Euler equations for the equivalent source method 
of viscous-inviscid interaction. Accounting for 
a permeable surface in the normal momentum rela- 
tion used to obtain pressure produced a slight 
improvement in the results. Finally, numerical 
solutions indicated some difference between the 
displacement surface method and the equivalent 
source method of viscous-inviscid interaction, 
although the difference in shock location was 
approximately the same as the distance between 
grid points. The eqivalent source method is the 
easiest to use and requires that only one grid 
be generated. 
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The present work was carried out on a CYBER 
203 In 6A-bit mode with only a small portion of 
the code vectorized. Typical run times were 208 
seconds for 1000 Euler equation cycles on a 
128 X 30 grid with 16 boundary-layer solutions. 
Experience indicates that a similar solution 
obtained on a CRAY-IS requires about half this 
amount of time. 
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